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£jTy Abstract. These pages review two families of mimetic finite difference methods: the 

^ | mixed-type methods presented in [3] and the nodal methods of [4]. The purpose of this 

exercise it to highlight the similitudes underlying the construction of the two families. The 
comparison prompts the definition of a piecewise linear postprocessing of the nodal mimetic 
finite difference solution, as it was done for the mixed- type method in [5]. 



1. Setting of the problem 

Let Q C M 3 be a Lipschitz bounded polyhedral domain. We consider the following elliptic 
boundary value problem: 

-V • (KVu) =g mn, 

u = in dfl, ^ ' ' 

where K G ~R dxd is a full symmetric tensor with components in W l,oc (Q). Strong ellipticity 
is assumed: thus there exists two positive constants and k* such that 

M v f < v T K(x)v < k*||v|| 2 Vv 6 R 3 , for a.e. xgO, (1.2) 

where || • || denotes the Euclidean norm of M 3 . 

Let Ll(fl) denote the space consisting of square Lebesgue-integrable functions having zero 
mean value and 

H{dw,Q) = {F e {L 2 {Q)) 3 : div F G L 2 {Q)} . 

In view of the discretization of fll.ip by both the nodal and mixed-type MFD, we consider 
the standard and mixed variational formulation of the problem, namely, 

find u e H]{n) : (KVu, Vv) = (g,v) Vi> G H^(n), (1.3) 

with (•, •) denting the L 2 -scalar product, and 

find (p,F) G L 2 (n) x H(div,n) s.t. 

(K- 1 F,G) + (p,divG) =0 VG G H(div,tt), (1.4) 
(divF,g) = (<7,g) Vg G L 2 (f2), 

respectively. 
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2. The mimetic finite difference method 

In this section we recall the mimetic finite difference methods introduced in [2], [3] and [I] . 
For more details we refer to the original papers and the references therein. 

Let Th be a sequence of non-overlapping conformal decomposition of Q into simply- 
connected polyhedral elements. For every element P we denote by \P\ its volume and by hp 
its diameter. Similarly, for each face / we denote by |/| its area and by hf its diameter, and 
for every edge e we denote by |e| its length. Depending on the context, dP denotes either 
the boundary of the element P or the union of the element faces. As usual, we set 

h = max hp. 

PeT h 

We assume the following mesh regularity conditions (see jl]). 

(HG) [Shape-regularity] There exist two positive real numbers N s and p s , independent of 
h, such that every mesh Th in the sequence admits a sub-partition Sh into tetrahedra 
such that: 

(HG1) every polyhedron P e % admits a decomposition Sh\p made of at most N s 
tetrahedra; 

(HG2) every tetrahedron T of Sh is shape-regular in the sense that the ratio between 
the radius rp of the inscribed sphere and its diameter hp is bounded from 
below by p s , i.e., 

^> Ps> WeS h . 
hp 

(ME) [Star-shaped elements] There exists a positive constant r*, independent of h, such 
that for each element P there exists a point Mp £ P such that P is star-shaped 
with respect to every point in the ball of center Mp and radius rjip. 

The MFD methods solution is a collection of real values associated to the set of elements, 
faces, edges, and nodes of the decomposition Th- Following [1], we thus introduce four 
discrete spaces V, T, S, and M made of collections of real values associated to each element, 
face, edge, and node, respectively. 

The value associated to a face or an edge is to be interpreted as flux or work of vector 
fields, and thus faces and edges are assumed to be given an orientation. For any face / we 
fix its orientation once and for all by attaching to it a normal vector rij. Further, any edge 
e with vertices (Vi, V 2 ), is assumed oriented from V\ to V 2 . 

The notation Mp will indicate the restriction of M to the nodes belonging to the element 
P, and so on. For any element P, we also introduce the notation Vp and fp for the number 
of vertices and faces, respectively. Further, we denote by V) the number of vertices of any 
face / 

The nodal MFD method discretizing (11 .3p is built over the set M of nodal values and 
the set S of edge values. The primal discrete differential operator is the gradient operator 
grad h : Af — > S defined as follows: for each Uh £ A/", given an edge e with vertices (Vi, V 2 ), 

(gr&d h u h ) e = ^-(u e 2 - u{), 
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where u\ denotes the value in Uh corresponding to the node Vi of the edge e. This is the 
discrete gradient operator introduced in [I], up to the scaling factor l/\e\: here we adopt such 
scaling for consistency with the definition of the divergence operator below (alternatively, 
we could have scaled instead the divergence operator, as in pQ). In order to implement the 
homogeneous Dirichlet boundary conditions, we shall also need to consider the subset Mq 
of the elements in M that are zero-valued on all vertices belonging to dVt. Similarly, the 
elements of the subset Sq of £ are zero-valued on all the boundary edges. Notice that the 
restriction to Ao of the operator grad^ maps into Sq. 

The mixed MFD method for solving (11.41) is built over the set V of elemental values and the 
set J 7 of face values. The primal discrete differential operator is in this case the divergence 
operator div h : J 7 — > V defined as follows: for each F h £ J 7 , given an element P, we set 

(div^)p= Yl l/l*f> 
' ' fedP 

where Fj = Ffiif ■ n P , with the normal of / out of P and Ff representing the value of 
F h associated to /. 

Let us remark that the two operators just defined operate on two different pairs of discrete 
spaces. In fact, it is possible to define, starting from the gradient and divergence operators 
above, a derived divergence and gradient operator, respectively. As these are not needed in 
the MFD formulation, we refrain to do so. 

We shall also need the relevant interpolation operators. Given any function p £ 
we define its interpolant H-pp £ V as 

(Jl P p)p = — [ pdV for all element P . (2.1) 
I "I Jp 

For every vector-valued function F £ (L S (Q)) 3 , s > 2, with divF £ L 2 (Q), we define its 
interpolant IT r/ £ J 7 as 

(Ujrf) f =l- fF-n f dS for all face/. (2.2) 
l/l J f 

As the interpolator for £ shall not be needed, we conclude by defining the nodal interpolator. 
Given any function u £ Hq(Q) fl C°(Q), we define its interpolant ILvu £ M as 

(ILvu)y = u(V) for all vertex V . (2.3) 

Similarly to finite element methods, an MFD method is defined by restricting a given vari- 
ational formulation to the MFD discrete set, with the crucial difficulty that the L 2 -product 
has to be substituted by discrete scalar products. Notice that, even when computing a finite 
element, the exact L 2 -product is substituted by a quadrature formula that is consistent with 
the degree of accuracy of the method. 

The principle used to define MFD methods is indeed that the discrete scalar product has 
to be consistent, i.e. exact on (the interpolants of) the correct space of polynomials: in the 
case of low order MFD methods, we shall require exactness on constants. Not surprisingly, 
the linear and Raviar-Thomas finite element methods are instances of the low order nodal 
and mixed-type MFD methods, respectively, obtained by a particular choice of MFD scalar 
products. 
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In order to derive the methods from the respective consistency conditions, we shall need 
to work with suitable approximations of the data. We denote by g the piecewise constant 
function obtained from the forcing function g by averaging over each element ? in 7/,. 
Similarly, IK will denote the tensor obtained from K by averaging each component over each 
P in %■ 

Up to the definition of the relevant scalar products (and linear functionals) , we write down 
the discrete MFD counterparts of the two problems ( II. 3p and fll.4p as follows. 
The nodal MFD method reads: 

find u h e Wo : [grad' 1 u h , grad h v h } £ = (g, v h ) M Vv h £ jV . (2.4) 

Further, we write the mixed MFD method as: 

find (F h ,p s h ) e FxV: 

[F h ,G h ] T +[p h ,div h G] v = VG ft GJ, (2-5) 
[div h F h , q h ] v = -[Uv9, Qh]v G ?. 

3. Scalar products 

The construction of the scalar products is achieved element by element and then summing 
up the elemental contributions. As mentioned earlier, the principle is that we want our 
scalar products to respect element by element the constants (more precisely, the interpolant 
of constant functions). 

The definition of the product in V is straightforward: 

Regarding the discrete space A/", we actually just need to define the linear functional 
(g, To this end, for each P e 7^, we introduce the numerical integration formula 



r Vp 

j vdp~j2 v ( v i 

j p i=i 



V P 

pjLUp, (3.2) 



where {wp}^ is a set of non-negative weights such that the quadrature is exact whenever 
v is a constant. We then define 

Vp 

(9, v h ) N ■= Yl 9\p Yl MV^oop. (3.3) 
PeT h i=l 

We now come to the less obvious problem of the definition of the scalar products mimicking 
the i^Q-products, namely [-, and [grad' 1 -,grad h -]s- 
The starting point is, in both cases, the Green identity 



$. V(j)dP= - / (f)V-^dP+ / <&-n P (j)dS (3.4) 
p Jp Jap 

valid for any sufficiently smooth vector-valued function $ and scalar function (f). 
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To deduce consistency conditions for the discrete scalar products, we specialize the above 
identity by testing it on the space of linear polynomials P 1 (P). To this end, let us consider, 
as basis of P 1 (-P), the set {bj}? =0 given by 

b (x) = 1, 

bj(x) = Xj ■ (x — xp), i = l,..., 3, 

where Xj is the j-th coordinate vector, and xp is the barycentre of P. 

Substituting </> = 6 m (E3J) gives back the divergence theorem, while with <fi = bj, j = 
1, . . . , 3, we obtain 

/ K- 1 ® ■ Kxj dP = - j bjV ■ <5>dP + j $ ■ n P bj dS, (3.5) 

J P JP JdP 

where we have also multiplied and divided by K in view of reproducing the weighted L 2 - 
product appearing in the first equation of (11.41) . 

This identity suggests to define, for every Gh G J 7 , and j — 1, ... ,3, 

[G^U^Kxj)}^ := f G P f b 3 dS 

f^ ap 1 (3.6) 

= E .no ■'■<■),(:';. 

f&dp 

with xj denoting the barycentre of the face /. Notice that the volume integral in the right- 
hand side of (13.51) disappears due to the fact that d\v h Gh is constant over P. 

It easily follows that any scalar product satisfying the so called local consistency condi- 
tion (I3.6P is exact on the interpolant of constant vectors. Indeed we get that 

[n^Kx^n^fcy^ = \P\%j = J^KVb, ■ Vb v (3.7) 

for all i, j = 1, . . . , 3. 

Substituting $ = KS/b in (|3T4|) gives the triviality = 0, while with $ = KVfy, j = 
1, . . . , 3, we obtain 

/ KVbj -V(j)dP = / KVbj ■ n P (j) ds = V (K P x, ■ nf) / <j>ds. (3.8) 
Jp JdP fedP Jf 

This time we mimic such identity at the discrete level by requiring that, for every Vh £ A/", 
and j — 1, ... ,3, 

[grad^K), grad^n^)^ = [v h , U^b^ := £ (K P nf ), ]T v{u>{. (3.9) 

fedP i=i 

Here, {wf}^ represents a set of non-negative weights of a quadrature formula used to 
approximate the integral over the face /. Assuming that such quadrature formula is exact 
for polynomials of degree < 1, we easily get that 

[U M b h ILjtfbjtff = \P\K hj = J KVh ■ Vfy Vi, j = 1, . . . , 3. (3.10) 
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Acceptable MFD scalar products over J 7 or M are bilinear forms satisfying respectively 
(13. 6p or (13. 9p which are symmetric and obey the following scaling properties: there exist two 
constants c* and c* independent of P G Th such that 

c ]T |P|(G^) 2 < [G/i, G h ]^ <c*J2 IPKGf) 2 VG h G .F, 

/e9P /eap 
c* £ ^(gradS^ < [v&jV&] p < c * \P\(&ad h v h ) 2 e Wv h G A/". 

eedP eedP 

The consistency, symmetry and positivity conditions leave some freedom in the definition of 
the forms, and indeed we shall get a family of MFD scalar products. This is better analyzed 
by considering the matrices associated to the bilinear forms. 

Let Mjr be the fp x f P symmetric matrix related to the form [-, ■]£ and the Vp x Vp 
symmetric matrix related to [•, ■] P f. 

We translate the conditions (13.61) and ( 13. 9 p into algebraic conditions for Mj? and Mfj by 
introducing the following matrices. Let N and R be the fp x 3 matrices given by 

N = [U^Kx,) . . . U r (Kx 3 )} = [n£ • • ■ nf fp ] T K P 

(the above equivalence is obtained expressing the definition of the interpolant Ujr) and 

\fi\(x fl -x P ) 

R= ; 

_ \fp\( x ff P ~ X P) _ 

respectively. We write the consistency condition ( 13. 6 j) as 

M^N = R. (3.11) 

Further, let W be the fp x Vp matrix collecting on each row the facial quadrature weights ap- 
pearing in ( 13 .9p filled with zeros to account for the nodes not belonging to the corresponding 
face. Then, by introducing the Vp x 3 matrices A and B given by 

A T = N T W and B = [ILv&i . . . ILv&a] , 

we write the consistency condition ( 13. 9 j) as 

iWj^B = A. (3.12) 

Let now C be a /y x (/y — 3) matrix with columns that span the null space of iV T and D 
be a Vp x (Vp — 4) matrix with columns that span the null space of [ILyfro ' ' ' Ha/-^]. Then 
the general form of acceptable matrices Mj% and Mfc is respectively given by 

M| = j^RK^RT + CU T C T and Af£ = -^AK^A 7 + DU N D T (3.13) 

with Ujr and Uj^ arbitrary symmetric and positive definite matrices of the appropriate scaling. 
The dimension of Ujr and Utf is (/y — 3) x (/y — 3) and (Vp — 4) x (Vp — 4), respectively. 
Notice that, if P is a tetrahedron, than we just get one possible nodal scalar product, namely 
M^f = jp-^AKp 1 A T . As noted in |4], when Th is made of tetrahedrons the nodal MFD values 

coincide with those of the standard P 1 finite element method (with K used in place of IK 
before evaluating the elemental integrals). 
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4. Gradient reconstructions 

It is clear that only the first part of the definitions of Mj and M^- given in (I3.13P acts on the 
relevant subspace of interpolated linear polynomials, and we know that such action is exact. 
As our forms operate at the gradient level, we can use them to define piecewise constant 
gradient reconstructions which have to be exact on linear polynomials. Following [5] , for any 
vector field G on P we define its reconstruction G R as 

Gf :=^[U T G,U T (Kx t )] R Vz = l,...,3, 

which is exact on constant fields by f !3.Tj) . Notice that the reconstruction is easily calculated 
as follows 

G R = j^N T M? (U^G) P = ^rK- 1 R T (U^G) p . (4.1) 

The reconstruction formula (14. ip was used in [5] to postprocess the mixed MFD method 
solution (p h ,F h ). Indeed, after (p h ,F h ) have been calculated, we can assemble elementwise 
a piecewise linear second-order accurate solution is given, on each P £ Th, by 

Ph \p = Ph\p + TpT (Kp^Ffclp) • (x - x P ). 

Similarly, we can define a reconstructed gradient of any scalar function v over P as 

(K P grad fl v)j := —- [Ilj^v , Uj^bj]^, Vj = 1, . . . , 3, 

which is exact on linear polynomials by (I3.10p . Also in this case we have the following 
formula for the computation of the reconstructed gradient: 

grad% = jLk-WmZ (U M v) p = y^^p 1 ^ (n^) P • (4.2) 

This formula may be used to define, starting from the nodal MFD solution u^, a piecewise 
linear solution u R given, on each P £ Th, by 

1 Vp 1 

u h\p=T^Yl u h(Vp)Up + 77^ (K p 1 A T U h \ P ^ ■ (X - X P ). 
' ' i=l II 

If the decomposition Th is made of tetrahedrons, the exactness of the elemental and facial 
quadrature formulas on the respective spaces of constant and linear polynomials implies that 
u^\p is the unique linear polynomial which takes the values Uh{V P ) at the elemental vertices. 
Thus coincides with the standard P 1 finite element solution. 

It is interesting to rediscover, in the case of tetrahedral elements, the exactness of (14. 2 j) on 
linear polynomials. To this purpose, let us consider {ifj}^ =l as the standard basis of P 1 (P), 
so that <fij(V P ) = Sij. Then, for any j = 1, ... ,4, denoting by fj the face opposed to the 
j-th vertex of the tetrahedron, we have 

grad^- = -gfpT l/il n £- 
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On the other hand, the reconstructed gradient of (fj is given by 

grad'V, = ^p 1 ^ (IW,-)p = ^(Kp 1 A T y = ±-([n^ ■ ■ ■ nfjWY, 

with the symbol (•) J ' indicating the j-th column of its argument. Now, the j-th column of W 
is given by (Wy = |(|/i| ■ ■ - • • ■ \ ff P \) T with the appearing in the j-th position, as these 
weights are the only ones ensuring exactness on linear polynomials. It follows that 

1 fp 1 1 
grad^ = —^-l/jlnj = -^|/>£. 

Pi 

5. Discretisation of advective terms: nodal MFD 
We now consider the advection-diffusion problem 

-V ■ (KVu) + -Vu = g infi, , . 

u = o in an, [b > 

where /3 is a given vector field with components in W '°°(Q). 

We want to construct a nodal MFD discretisation of (15. ip . The discretisation of the new, 
advective, term will be based on the gradient reconstruction formula (14. 2 j) and on a piecewise 
constant approximation of the data. 

Let P be the vector field obtained from /3 by averaging each component over each mesh 
element P. Then /3 ■ grad h % e V, for all v\ t £ Af. We define the following nodal MFD 
method: 

find u h G A/o : [grad h grad h u/Jg + (/3 • grad R u h , v h )x = (g, v h )x V v h e A/" . (5.2) 

We may as well introduce in the method the following streamline-diffusion type stabilizing 
term: 

[f h j3 ■ grad R u h , j3 ■ gmd R v h \ v , (5.3) 

where fh is a stabilization parameter defined element-wise in function of the local mesh 
Peclet number. 

Notice that, on each element P, the reconstructed gradient appearing above gets the 
following expression: 

grad^^lp = j p - l Kp 1 A T v h \ P 

11 Sp y (5.4) 

1 1 fedP i=i 

Thanks to the exactness on linears of the gradient reconstruction, the term (I5.3p satisfies 
a local consistency condition similar to (13. 9p . Once again, the starting point is a Green 
identity. Given a generic smooth function and a linear function p G P 1 (P), we have 

/ • Vp)(/3 • V<P)dP = [ ■ Vp)09 ■ n P )^dS = ■ Vp) V (/3 • n/ ) / (5.5) 
Jp Jap /eSP Jf 



The term (15.31) mimics ( 15. 5 P in that, for each Vh G M and p G P 1 (P), 



\P ■ grad fi IVp, f> ■ grad R = \P\0 ■ grad fl ILvp) (0 • grad fl v h ) 

Vf 



(5.6) 



c/. Equation (5.14) in Indeed (15. 6p follows from (15. 4p and the fact that grad ITa/p = Vp, 
and thus (3 ■ grad R Hj^p = (3 ■ Vp. Actually, this last equality can be re-obtained by direct 
calculation: again from ( 15 .4p we have that 
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